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Abstract 

The modification of Smoothed Particle Hydrodynamics (SPH) method with Riemann 
Solver is called Godunov SPH. We further extend the Godunov SPH to the description 
of a medium with negative pressure. Under certain circumstances, the SPH method 
shows an unphysical instability that results in particle clustering. This instability is 
called the tensile instability. The tensile instability occurs in positive pressure regions in 
a regular fluid if a very large number of neighbor particles are used with certain shapes of 
kernel functions, and it is significant in negative pressure regions that emerge in stretched 
elastic bodies. We must suppress the tensile instability in SPH for calculations of elastic 
bodies. In this study, we develop a new technique to remove the tensile instability by 
extending the Godunov SPH method and conducting a linear stability analysis of the 
equation of motion for the extended method. We find that the tensile instability can be 
suppressed by choosing an appropriate order of interpolation in the equation of motion 
of the Godunov SPH method. We also derive an analytic solution for a Riemann solver 
for a simple equation of state of an elastic body, and construct a Godunov SPH method 
for the equation of state that allows negative pressure. 

Keywords: Smoothed Particle Hydrodynamics, Tensile instability. Linear stability 
analysis, Godunov’s method 


1. Introduction 

Smoothed Particle Hydrodynamics (SPH) is a computational fluid dynamics method. 
It is a Lagrangian method and does not require a Eulerian mesh. Each SPH particle 
mimics a fluid element and we can describe fluid dynamics by the motion of SPH particles. 
Lagrangian particle methods like SPH are suited for systems that have large voids or large 
deformed structures. Thus, the SPH method has been widely used in astrophysics and 
planetary science since it was proposed by Lucy [H and Gingold & Monaghan Q. The 
most popular form of SPH is called the standard SPH method. 

Recently, some studies have attempted to apply SPH to elastic dynamics (e.g., ii 
H). The application is straightforward because the equations of elastic dynamics are very 
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similar to those of hydrodynamics. Furthermore, we can easily follow the crack history 
of elastic body elements because SPH is a Lagrangian particle method. Therefore SPH 
is a powerful tool for the calculation of disruptive collisions. 

However, the standard SPH method for elastic dynamics has one serious problem. 
As opposed to a compressed elastic body, which produces positive pressure, a stretched 
elastic body produces negative pressure. In the negative pressure regime, the standard 
SPH method has an instability that results in a clustering of particles. This instability 
is called the tensile instability, and was studied in detail by Swegle et al. It is 

also known that this tensile instability occurs even in the positive pressure regime. Q 
demonstrates that B-spline kernels produce the tensile instability, even in the positive 
pressure regime, if the number of neighbor particles is large. If the number of neighbor 
particles is small, we can dodge the tensile instability in the positive pressure, but this 
increases the ‘i?o error’ [^, which is the leading error of momentum equation, and will 
make erroneous result. Other discretizations such as Q successfully reduce this ‘i?o 
error’, but these discretization forms lose momentum conservation. These facts may 
suggest that the conservation and reordering property are not easily compatible with the 
method that does not suffer from the tensile instability. This instability is unphysical. 
Thus, for example, when we calculate collisions of asteroids without any prescription for 
the tensile instability, the size distribution of fragments may become completely wrong. 
Therefore, it is important to develop a method that suppresses the tensile instability. 

There are many approaches toward the solution of the tensile instability in the stan¬ 
dard SPH method (e.g., 13 ) 11 > 12 ; 13 ) 3 13 )• In 121, Monaghan introduced an artificial 
pressure that provides a strong repulsive force only when particles become close to each 
other. He also conducted a linear stability analysis for this method and found that this 
artificial pressure can suppress the instability at short wavelengths but does not affect 
the perturbations of long wavelengths. However, according to Mehra et al. 1^, this 


artificial pressure cannot suppress the tensile instability in simulations of hypervelocity 
impacts. 

Another formulation of SPH is called the Godunov SPH method [ 13 . It achieves the 
second-order spatial accuracy, whereas the standard SPH method does not have such a 
convergence property because of its rough approximation. The tensile instability does 
not occur if we solve the original equations of hydrodynamics or elastic dynamics exactly. 
Thus, this instability is caused by discretization error. We expect that if we use a method 
that has higher-order accuracy, such as the Godunov SPH method, we can suppress the 
tensile instability. 

In this paper, we formulate a higher-order interpolation of the equations of the Go¬ 
dunov SPH method. Furthermore, we evaluate the stability of each order of interpolation 
by a linear stability analysis. A practical approach for dealing with the tensile instability 
using variable smoothing lengths is also shown. We also derive the analytical solution 
of the Riemann problem for a simple equation of state for an elastic body that allows 
negative pressure, and we construct a Godunov SPH method for this equation of state. 

The structure of this paper is as follows: in Section 2, the essence of the SPH method 
and the equations for the Godunov SPH method are introduced. In Section 3, we present 
the motivation for higher-order interpolation, and then we derive the equations for higher- 
order interpolation. We also derive the analytical solution of the Riemann problem 
for a simple equation of state for an elastic body. In Section 4, we conduct a linear 
stability analysis of the equation of motion for the Godunov SPH method and evaluate 
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its effectiveness for the tensile instability. Test calculations are presented in Section 5 
and we show the validity of the analysis presented in Section 4. Section 6 is a summary 
of our work. 


2. SPH Method 


2.1. Essence of SPH 

The contents of section 2.1 and 2.2 follow section 2 in 
Lagrangian particle method, and hence we use the Lagrangian forms of the equation of 
motion and equation of energy for an inviscid fluid: 


17|. The SPH method is a 


dv 

dt 

du 

dt 


--yP, 

P 

P^ 

-V • V, 

P 


( 1 ) 

( 2 ) 


where u is the specific internal energy, v is the velocity, P is the pressure, and p is the 
density. 

In the SPH method, a physical quantity, /, at an arbitrary position is approximated 
by the convolution of nearby quantities. The convolution of a quantity / at position r 
is defined as 


{f){r) = J f{r')W{r'-r,h)dr', (3) 

where W{r — r,h) is a kernel function and /i is a parameter called the smoothing 
length. We use the angle brackets () to denote the convolution. We tentatively treat this 
smoothing length as constant. 

The kernel function has various forms; one of the simplest kernels is Gaussian, 


W{r, h) 



(4) 


where d represents the spatial dimension. We use this Gaussian kernel throughout this 
paper. 

Now we define the density at arbitrary position by the summation of the kernel 
function at particle positions, 


Pir) = J2mjW{r-rj,h). (5) 

3 

From this we can make the following identity: 


i = E 

3 


^3 

p{r) 


W{r 


rj,h). 


( 6 ) 


Using the identity given by Eq. dH]), and the definition of the convolution, Eq. ([31), we 
can express a physical quantity at particle position i as 
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( 7 ) 


fi={f)in) = J fir')W{r' -n,h)dr' 

= f 'mj^-j^W{r — ri,h)W{r —rj,h)dr . 

Similarly, we can express the space derivative of a physical quantity of particle i as 

^ (fr “ n,h)dr' 

= Zl/ -r,,h)W{r' -rj,h)dr', (8) 


where 


d / d d d \ 

dvi V dxi ’ dyi ’ dzi J' 


(9) 


2.2. Standard SPH 

In the previous subsection we introduced the general formalism of the SPH method. 
However, the formalism that is used in the standard SPH method is simpler than Eq. © 
and Eq. ([5]). 

In the standard SPH method, we integrate Eq. © using the approximation IE(r — 
rj,h) Ki 6{r — rj), such that 

fi!=z^mj—W{r,-rj,h), ( 10 ) 

where 5{r) is the Dirac <5 function. Similarly, we can integrate Eq. ([8|) and express the 
space derivative of a physical quantity of particle i as 

j \. J I 

Pj i 
3 

These expressions are derived from rough approximation, but in actual calculations 
Eq. (fTTI) is sufficient for the spatial derivative of a physical quantity. 

2.3. Godunov SPH 

In this subsection, we introduce the equation of motion and the equation of energy 
that were derived by [l7| . We also introduce the equations for variable smoothing length 
in the Godunov SPH method. 
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2.3.1. Equations for Godunov SPH 

The equation of motion for the Godunov SPH method is defined by the convolution 


of Eq. (P) 



( 12 ) 


where the overdot shows a time derivative. 

We transform the right-hand side of Eq. (fT^ using Eq. ([5]) and integration by parts, 



(13) 


If we multiply both sides of Eq. (1131) by the mass of particle i, the left-hand side 
becomes the time derivative of the momentum of particle i, and the right-hand side 
becomes antisymmetric with respect to i and j. Therefore, in this formalism, as in 
standard SPH, the linear momentum and angular momentum of a particle system are 
conserved. 

When the equation of state only depends on density, or the fluid is barotropic, we 
can follow the evolution of the fluid using only the equation of motion. However, in the 
case of a general equation of state like P = P{p,u), or a fluid that has a shock wave, 
the equation of energy is also required. The equation of energy for the Godunov SPH 
method is derived from the convolution of Eq. (P, similarly to the derivation of Eq. (Ha, 



(14) 


We integrate the right-hand side of Eq. m by parts and find 



Then, we use the following approximation: 



(16) 


By using Eqs. ([TC| and (ini) we can transform Eq. (IlSp into 



Equations o and (dzi) are not yet useful for practical calculation because they 


contain spatial integration. Thus further approximations are required. 
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2.3.2. Convolution 

The equation of motion and the equation of energy for the Godunov SPH method 
are shown in Eq. m and Eq. CZD, and these equations involve spatial integration. The 
p{r) parts of the integrands are given by Eq. ([5]). Thus, we cannot integrate Eq. (fT51) and 
Eq. (1171) analytically. Furthermore, it is almost impossible to integrate these equations 
numerically for each pair of i-th and j-th particles in practical calculations because of the 
heavy numerical cost. Therefore, we must interpolate and assume the spatial distribution 
of physical quantities like p{r). 

We now define a new convenient coordinate system for the integration. This coor¬ 
dinate system has its origin at (r^ -|- rj)j2, and we define the s-axis to be along the 
direction of the vector — Vj. We use sj_ to denote the component of vector r that is 
perpendicular to the s-axis, and we define e^- = (r^ — rj)/\ri — rj\ as the unit vector 
along the s-axis, and As^ = — rj\ as the distance between particles i and j. 

Essentially, Eq. (1131) and Eq. (1171) include the integration as follows: 

/ - ri,h)W{r - rj,h)dr. (18) 

J P^[r) 

For simplicity we define the weighted average /j* as 


[ Vi, h)W{r-rj, h)dr = f*j f ^—W{r - ri,h)W(r - rj,h)dr. (19) 

J P\r) J p^{r) 

Then we expand p~^(r) linearly in the direction perpendicular to the s-axis as 

p“^(r) Ri p“^(s)-f s_L • Vp“^(s). (20) 

Note that if we use the approximation of Eq. (EHl) the component that is perpendicular to 
the s-axis vanishes because of the symmetric property of the kernel function. By using 
Eq. (Oni) we can transform Eq. dT51) into 


where 


j “ ri,h)W{r - rj,h)dr = - r^, v^h), 

1 




p^{s) 


exp 


(-f) 


ds. 


( 21 ) 


( 22 ) 


Finally, if we interpolate p{s) in the direction of the s-axis, we can integrate 17^ 
analytically. Here it is convenient to define the specific volume and its gradient as 


V{r) = 
VV{r) 


1 

p(r)’ 


1 

p2(r) 


Vp(r) = 


1 

p2(r) 


^mjVW(r-rj,h). 

j 


(23) 
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The simplest interpolation of l^(s) is linear interpolation, which we express as 




A 


(24) 


where 


= 


Dij — 


Vjn) - Vjr,) 

Asij 

V{ri)+V{rj) 


(25) 


By substituting Eq. into Eq. (12^ . we obtain 

1 


V^,,nneM = -h^Cf^+Df^. (26) 

Another convenient interpolation is cubic spline interpolation. This is expressed as 


E(s) — —y — AijS^ + BijS^ + CijS + Dij, 


(27) 


where 


, E, +E- 


A - ■ — —9 

“ AsT ' As2. 

y y 

1 V- - V- 

TD _ ‘ J 

" 2 As,, ’ 


A, = + E,) - i(E; - e;)As,„ 


E = E(r,), 

E- = E(r,), 

E; = e,, • VE(n), 
E; = e,, • VE(r,). 

Then we can calculate Eq. (l22ll and obtain 


(28) 


E,^,cubic(/i) = + 52) + \h^mjD,3 + C^) + 5^. (29) 

If we substitute Eg. (051) or Eg. (1^ and Eg. (1^ . into Eg. (fOl) and Eg. (fT71) . the 
equation of motion and the equation of energy finally become 
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Vr = -2j2mjP*jV^jih) — W{ri - rj,V2h), 


( 30 ) 


= - P:^n) ■ V?,(h)^W{r, - r„V2h). 


(31) 


In 17[, Inutsuka introduced a Riemann solver that uses the physical quantities of 
the i-th and j-th particles as initial conditions of the shock tube problem. We use this 
Riemann solver for P*^ and v*j and can introduce necessary but minimal viscosity to 
deal with shock waves. However, in this paper, in order to separate the effect of viscosity 
and the effectiveness of the formalism of the Godunov SPH method against the tensile 
instability, we use Eq. (1!?^ for P*j when we analyze the stability. 


P* = 
^ij 


P,. 


Pj 


(32) 


If we use Eq. (15^ instead of the Riemann solver, we can turn off the viscosity. 

Throughout this paper we use the following simple equation of state (e.g., [H): 


P = Cf(p-po.eos), (33) 

Thus, we do not use the equation of energy given by Eq. dlTTI) . Here, po^eos is a reference 
density and Cs is the sound velocity, which we treat as a constant. 


2.3.3. Variable smoothing length 

We have so far treated the smoothing length as constant. The smoothing length 
should be close to the average particle spacing, and the average particle spacing varies 
largely in space when the density varies largely. In practical calculations we should vary 
the smoothing length according to the local mean of the particle spacing. In this section, 
the smoothing length is represented by spatial variable h{r). Then the integration of 
Eq. (IT!?!) and Eq. (ITTll includes /i(r), and we cannot integrate analytically even if we use 
the polynomial approximation of p~^{r). 

In |l7l| , Inutsuka used approximate analytical integration assuming that the smooth¬ 
ing length is hi for the half of the integration space that includes particle f, and hj for 
the other half. Then the equation of motion and the equation of energy for variable 
smoothing length become 


Vi=-Y^ mjP* - rj,V2hi) 


= - Ptj'^i) ■ - rj,V2hi) 


(34) 


+ ^^“ ’’I ’ ) 


(35) 
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This approximation assumes that the smoothing length does not vary much within the 
neighborhood of particles i and j. 

A possible way to decide the smoothing length of particle i is, for example, 


where 



(36) 


p* = '^mjW{ri - Tj.h*), h* = hiCsmooth, (37) 

i 

and Csinooth is a constant to smooth the spatial variation of the smoothing length, p 
is a constant to determine the ratio of the smoothing length to the average particle 
spacing. In this study we use 77 = 1. Csmooth > 1 we can obtain smoother variation of 
the smoothing length than that of the density, and the approximation that was used to 
derive Eq. (13411 and Eq. (1351) becomes more accurate. However, the number of neighbor 
particles to obtain correct smoothing length becomes large if Csmooth becomes large. 
Eq. (1551) and Eq. (1571) are recursive equations that require iterative calculation of hi . In 
practical calculations the use of hi in the previous time step for hi in Eq. (13711 works well. 


3. Further Extension of Godunov SPH 

In this section, we show that the original formalism of the Godunov SPH method is 
in principle stable against the tensile instability. Next, we derive an equation for higher- 
order interpolation of 14^ for use in Eq. (1151) . We also derive an analytical solution of 
the Riemann solver for Eq. (1551) so that we can use the Godunov SPH method for an 
equation of state that allows negative pressure. 

3.1. Motivation 

The equation of motion for the Godunov SPH method is given by Eq. (1151) before we 
interpolate the distribution of physical quantities. To derive this equation we only take 
the convolution of the original equation of motion for an inviscid fluid. This convolution 
only introduces smoothing of the physical quantity, and thus should not cause instability. 
This suggests that the exact integration of Eq. o should remove the tensile instability 
even in the case of negative pressure. Note that the density at arbitrary position is given 
by Eq. ([S]) and the pressure at arbitrary position is calculated from Eq. (1551) . Thus we 
can calculate the convolution integral in Eqs. (1131) and (1221) using numerical integration. 

To show this, we derive the dispersion relation for Eq. m- For simplicity, we focus 
on the one-dimensional case and calculate the dispersion relation almost exactly using 
the method shown in Appendix A. In short, we assign particles with equal spacing and 
add a sinusoidal perturbation with wave number k and calculate the acceleration. Then 
we can calculate the squared angular frequency by taking the ratio of the displacement 
and acceleration. 

Figure [T] shows the dispersion relation of Eq. (IT51) that is calculated using the method 
of Appendix A. Here we consider all particles to have equal mass, and the parameters are 
as follows: mass of each particle m = 0.0008, particle spacing in the unperturbed state 
Ax = 0.001, h = 0.001, Cs = 1.0, po.eos = 1-0. The average density is po = 0.8 and the 
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average pressure is Pq = —0.2 (negative). Numerical integration is done using a simple 
trapezoidal law, and the integral interval is Ax/10. The unit of length is normalized by 
the average particle spacing. Thus, k = tt corresponds to a perturbation wavelength of 
two particle spacings. Negative means the perturbation of the given wave number is 
unstable. 



k 


Figure 1: Dispersion relation of Eq. l|T3ll for negative pressure. The horizontal axis is k, the vertical axis 
is and the unit of length is normalized by the average particle spacing. Solid curve shows the 

exact dispersion relation of the sound wave. 


As shown in Fig.[Tl > 0 for all wave numbers. Thus, we can see that the exact 
integration of Eq. m can remove the tensile instability for negative pressure. 

Next, we examine a simpler case. In this case we use the equation of motion given 
by Eq. (1501) . but we integrate Eq. (I00|) numerically. Here we use P*j given by Eq. (I50|) . 
This is the same as when we integrate Eq. (USD, but the pressure is assumed constant, 
P{r) = {Pi + Pj)/2. Eigure [0] shows the dispersion relation when we calculate the 
acceleration using this method, where all parameters are the same as in the case of 
Fig.ffl 

In Fig.jO] it appears that = 0 at around fc = tt, but actually these are small 
value and not 0. Thus, our calculations are free from the instability, even for negative 
pressure, if we integrate numerically. 

As shown in the test calculations in section 5, the Godunov SPH method also suffers 
from the tensile instability in the case of inappropriate interpolation methods. Therefore, 
the exact integration of the convolution can suppress the tensile instability. These results 
imply that the tensile instability of the Godunov SPH method comes from the errors in 
the approximation for integration. Thus we expect that if we improve the interpolation 
of Fjj we can remove the tensile instability. 

3.2. Quintic spline interpolation 

As we discussed, if we use a higher-order interpolation of to achieve closer to 
the exact integration of Eq. (1551) . it is possible to suppress the tensile instability. In this 
subsection, we derive an equation for the higher-order interpolation. 
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Figure 2: Dispersion relation of Eq. m for negative pressure, where Eq. is integrated numerically. 


In [17[ , Inutsuka constructed a cubic spline interpolation of V(s) by using four quan¬ 


tities: the specific volume and its derivative for particles i and j. We can construct a 
quintic spline interpolation by adding two more quantities: the second derivative of the 
specific volume for particles i and j. 

To derive the second derivative of the specific volume for an arbitrary direction, we 
simply use Eq. m and differentiate the first derivative to obtain 


d^Vjr,) 

dxadxp 


E 


rrij dV (rj) d 

Pj dXa^j dxp^i 


W[ri 


rj,h), 


(38) 


where a and /3 take values 1, 2, or 3, and xi = x, X 2 = y, x^ = z. Actually, the 
second derivative that is calculated with Eq. (l38l) is not so accurate and we obtain a 
smoother second derivative than the exact value, but this method is sufficient for the 
tensile instability as we will show using the linear stability analysis of Section 4. In the 
actual simulations, we calculate the second derivative using Eq. (1381) after calculating the 
first derivative using Eq. (1^ . 

Once we calculate the second derivative, we can construct the quintic spline interpo¬ 
lation of V{s). This is expressed as 


V(s) — AijS^ BijS^ CijS^ DijS^ EijS Fij^ (39) 

where 


II 









and 




- Vj _ 3 .^* +^j , 


+ 


As'^- ' 2 As^- 

y y 


IV, -V^ , ly, 

-t>7.'>- = — - 


2 Aa^ 


4 A4 


.Vi-V 5 Id + K- 1 Vi - V 


^'Asl ' 2 As?, 4 As, 


3 1 / " " \ 

^« = iA;r■ s*’"- -"d’' 

= f ^ +d >+^(V - d>»«. 

a, = ^(''i + Vj) - ^(d - vdAsi, + ^(d' + d')A»«. 


(40) 


^ ®y.«®y,/3 

a p 

^j ~ 6y.a6y,/3 

a P 


d^Vjr,) 

dXadxp 

9^V{rj) 

dXadxp 


(41) 


Here, eij^a shows the a direction component of . By substituting Eq. (l3^ into Eq. 
we hnd 


V' 


2 j,quintic 


) + ^h<^(2A,,E,, + 2B,,D,, + Cf^) 
+ ^h\2B,,F,, + 2C,,E,, + Dl) + \h\2D,,F,, + E^) + . (42) 


3.3. Riemann solver for elastic equation of state 

Realistic equations of state generate negative pressure under certain circumstances. 
To develop a Godunov SPH method for such an equation of state we should derive the 
appropriate solution of a Riemann problem that includes shock waves with negative 
pressure. The emergence of a shock wave with negative pressure may not be obvious. 
As explained in many textbooks (e.g., 18|), a shock wave is generated by nonlinear 
steepening of a sound wave of finite amplitude. This nonlinear effect exists even in 
negative pressure media, and we have confirmed steepening in a negative pressure state 
for Eq. (Fig. [ 13 ] of Section 5). This fact implies the emergence of a shock wave with 

negative pressure. Therefore, it is valuable to construct a Riemann solver for such an 
equation of state and develop a Godunov SPH method that can handle shock waves in 
negative pressure regions. 
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In this subsection, we derive the analytical solution of a Riemann problem that uses 
the equation of state of Eq. (I33|). A method to solve a Riemann problem is called a 
Riemann solver. The Riemann solver for an ideal gas has been extensively analyzed 
in 0- Godunov SPH uses the resultant pressure P* and velocity v* of the Riemann 
problem for P*^ and v*^ of the equation of motion, Eq. (HHH) , and the equation of energy, 
Eq. (|3I|). 

Let us consider a shock wave propagating with velocity Vs in a medium of density 
p moving with velocity v. According to the Rankine-Hugoniot relations, we have the 
following relations across the shock wave: 


±Ws{y*-V) + {v*-v) = 0, (43) 

±Ws{v*-v) - [P* - P) = {), (44) 

where Ws denotes the Lagrangian shock speed |/9Us|, the sign shows the direction of 
shock propagation, and the quantities marked by asterisks correspond to the post-shock 
values. Eq. (14311 and Eq. (1441) represent conservation of mass and momentum flux across 
the shock wave. Using Eq. (H5|) . Eq. (H^ . and the equation of state for the post-shock 
value. 


P* = - P0,eos), 

we can express Ws as 


(45) 


Ws=Cs 



P* - P\ 

~c^)- 


(46) 


Next, to derive the relation between the physical quantities across the rarefaction 
wave, we derive the Riemann invariant for the equation of state given by Eq. (1331) . The 
Eulerian continuity equation and equation of motion in one dimension are 


dp dv dp 

dv dv 1 dP 

dt~^^dx p dx' 

Using Eq. (IHHll . we can transform Eq. (HHll into 

dv dv ^ dCs In p 

If we multiply both sides of Eq. (ITTl) by Cslp, we obtain 

dCs In p . ^ dv dCs In p _ 

Finally, by taking the sum and difference of Eq. (H^ and Eq. (1501) . we obtain 
— {v±Cslnp) + {v±Cs)-^{v±Cslnp) = 0 . 

13 ^ 


(47) 

(48) 

(49) 

(50) 


(51) 







Equation (EB means the Riemann invariants J± = ± Cg In p are conserved on 

a trajectory dx/dt = u ± Cg. These Riemann invariants remain constant across the 
rarefaction wave, thus 


r; =F Cg In p = n* =F Cg In p*, (52) 

where the sign shows propagation in the opposite direction to the rarefaction wave. If 
we define Wg for the rarefaction wave as in Eq. (H^ . and use Eq. (l5^ and Eq. (HSl) . we 
obtain 


ITg = 


p* -p 

P-P* 

hn( 

-1 

— V 

Cs 

V P* -|- Gg Po,eos 



(53) 


We can calculate P* and v* using Eq. (l46ll and Eq. (l53ll . Erom the definition of Wg, we 
obtain the following relations for the waves that propagate in the right and left directions: 


-Ws^l{v*-vl)-{P*-Pl) = 0, (54) 

- vn) - {P* - Pn) = 0, (55) 

where L and R denote the values on the left-hand side and right-hand side of the initial 
discontinuity, and 


Ws,B = Cg^PD^|OD-l- ^2 


VEg^D = 


Pn-P* 


Cs 


In 


C^PB 


P* + C^Po, 


if P* < Pd, 


(56) 


where D represents L or R. Eliminating v* or P* from Eq. (IMl) and Eq. (I5S|) , we obtain 


,* _ PhlWs^L + Pk/Ws^K + 

“ l/VEg,L + l/VEg,R 

H _ VhWs^L + fRlTg^R + Pp — Pr 

IEg,L + bbg.R 


(57) 

(58) 


We can determine P* and v* by iterative calculation of Eq. (1551) and Eq. (1571) . In 
practical calculations, 5 cycles of iterations with the initial value P* = {Pl -I- Pr)/2 are 
sufficient. 

Many general equations of state, such as the Tillotson equation of state (e.g., @), 
are composed of a combination of an ideal gas equation of state, a polytropic relation 
P = Kp'*, and the equation of state given by Eq. (1551) . Therefore we expect that a 
Riemann solver for general equations of state can be written as a combination of these 
Riemann solvers. 


4. Result of Linear Stability Analysis 

In this section, we conduct a linear stability analysis of the Godunov SPH method 
and evaluate its stability against the tensile instability for linear interpolation, cubic 
spline interpolation, and quintic spline interpolation of 14^. 
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4..1. One-dimensional case 

First, we conduct a linear stability analysis of sound wave propagation using a one¬ 
dimensional code. Figure [3] shows the dispersion relation for linear interpolation, cubic 
spline interpolation, and quintic spline interpolation in a one-dimensional flow with nega¬ 
tive pressure. The parameters are the same as in Fig.[T]and Fig. [31 m = 0.0008, h = 0.001, 
Aa: = 0.001, Cg = 1.0, po.eos = 1-0, and the average pressure is Pq = —0.2 (negative). 
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Figure 3: Dispersion relation for a sound wave with negative pressure using linear interpolation (red 
squares), cubic spline interpolation (green circles) and quintic spline interpolation (blue triangles) in a 
one-dimensional code. The points show the results obtained using the method shown in Appendix A, 
and the solid curves show the analytical solution for each interpolation. 


Each point shows the results that are calculated using the method shown in Appendix 
A, and the solid curves show the analytical solutions that are obtained from a linear 
analysis of each interpolation. Here, we only show the formula of the analytical solution, 
and the details of the linear analysis are shown in Appendix B. 

In the case of linear interpolation. 


where 


UJ 


2 

linear 


-C^gDa + 


2DPo ^ 

- Cl + 

Po 


2Po, 

- 0 . 

Po 


(59) 
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In the case of cubic spline interpolation, 


where 
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Po 


d = ^^{xi — Xj){l — cos[fc(a;i — Xj)])-gzz'^ixi — Xj,V2h) —. 
In the case of quintic spline interpolation, 


(62) 
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where 


An = '^{xi - Xj)" sin[fc(a;i - Xj)\-^W{xi - Xj,V2h)^, 




Po 


Bn = '^{Xi - a;j)”(l - cos[A:(a;i - Xj)])-^W{xi - xj, V2h)^, 




Po 




C = y^(l - cos[fc(a;i - Xj)])-^^W(xi - Xj, \/2h) —, 


(64) 


and Uci and UcJ show the position of a particle in the unperturbed state. 

As we can see in Fig. [31 the dispersion relations obtained using the method of Ap¬ 
pendix A agree with the analytical solutions given by Eq. (1531) , Eq. (1511) , and Eq. (153)) . 
In the cases of linear interpolation and cubic spline interpolation, < 0 at short wave¬ 
lengths (large k). Thus these calculations are unstable. However, in the case of quintic 
spline interpolation, > 0 even at short wavelengths, and so this quintic spline inter¬ 
polation is stable for negative pressure. 

These results do not depend on the absolute value of Pq or po, as long as the pressure 
is negative, because a in the first term of Eq. (1531) and Eq. (1511) , and Aq in the first term 
of Eq. (1531) . become 0 at the Nyquist frequency (fc = 27r/2Aa;) and all the other terms are 
proportional to Po/po- Therefore, these results might be general as long as the pressure 
is negative in the one-dimensional case. 

Surprisingly, this fact shows that, at least at the Nyquist frequency, > 0 for 
negative pressure implies < 0 for positive pressure. Figure |3| shows the dispersion 
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relation for the same conditions as in Fig. [3] but with po,eos = 0.6 and hence positive 
pressure, Pq = 0.2. 
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Figure 4: Same as Fig. |3] but for positive pressure. 


As shown in Fig.|4l linear interpolation and cubic spline interpolation are stable, but 
quintic spline interpolation is unstable at short wavelengths. Therefore, we conclude 
that quintic spline interpolation is stable for negative pressure but unstable for positive 
pressure in the one-dimensional case. 

4-2. Two-dimensional and three-dimensional case 

Next, we conduct a linear stability analysis for the two-dimensional and three-dimensional 
cases. We put the particles on a Cartesian lattice with a side length of Ax in the unper¬ 
turbed state. For simplicity we assume that the wave number vector is along the x-axis. 

In this case the analytical solutions of the dispersion relation for linear interpolation and 
cubic spline interpolation are 



(65) 


( 66 ) 


Actually, the forms of Eq. (1551) and Eq. (1551) appear to be the same as the one¬ 
dimensional case, Eq. (1551) and Eq. (15lT) , but the definitions of the coefficients are different 
and include summation for the y- and z-direction; 
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dxi 
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Po 


d = ^^{xi — Xj){l — cos[k{xi — Xj)])-gZzW{ri — rj,\/ 2 h) —. 
j^i ® do 


(67) 


Note that the values of the coefficients, apart from c, are almost independent of the 
number of the spatial dimension. The dispersion relation for linear interpolation does not 
have coefficient c, and the results of the two-dimensional and three-dimensional cases are 
the same as those of the one-dimensional case; linear interpolation is stable for positive 
pressure and unstable for negative pressure. 

The analytical dispersion relation for quintic spline interpolation in the two-dimensional 
and three-dimensional cases are 


<^quintic — ~ C^DAo.o H C — 


Po 


— \-\h^CpBi^^ - ^h‘^CpDA2,4 + lh^CpB4,2 
Po L 8 16 4 
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( 68 ) 


where 


An,m = £_^|L sin[/c(a;i - x,)\-^W{r, - v^h)^, 


Po 




Po 


92 


C = ^(1 — cos[fc(a;i — Xj)])—^W{ri — rj,V 2 h) 

OXi 


Po 


(69) 


Figure [ 5 ] shows the dispersion relations for cubic spline interpolation and quintic spline 
interpolation in the two-dimensional and three-dimensional cases for positive pressure, 
and Fig.ini shows the same dispersion relations but for negative pressure. Here, the 
parameters are the same as the one-dimensional case (Fig. [3] and Fig.H]), but to set the 
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Figure 5: Dispersion relation for cubic spline interpolation in the two-dimensional case (red squares), 
cubic spline interpolation in the three-dimensional case (green circles), quintic spline interpolation in 
the two-dimensional case (blue triangles), and quintic spline interpolation in the three-dimensional case 
(black diamonds) for positive pressure. 
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Figure 6: Same as Fig.^but for negative pressure. 


average density to po = 0.8 we set the mass of each particle to m = 0.8 x 10“® in the 
two-dimensional case and m = 0.8 x 10“® in the three-dimensional case. 

As we can see in Fig.[5]and Fig.|6l quintic spline interpolation has the same result as in 
the one-dimensional case; it is stable for negative pressure and unstable for positive pres¬ 
sure. However, the correspondence is very different for cubic spline interpolation. Cubic 
spline interpolation has the opposite result to the one-dimensional case; it is unstable for 
positive pressure and stable for negative pressure. 

This difference between the one-dimensional case and the two- and three-dimensional 
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cases for cubic spline interpolation comes from the coefficient c, because all the coeffi¬ 
cients, except for c, are almost the same for all spatial dimensions. 

4.3. Variable smoothing length case 

The extension of the method to the variable smoothing length case is not straight¬ 
forward because Eq. (IMl) is not exact in this case. However, we can satisfy almost the 
same condition as for constant smoothing length in the neighborhood of the f-th and j-th 
particles if Csmooth becomes sufficiently large, because it produces very smooth spatial 
variation of the smoothing length. Thus we expect that if Csmooth becomes large we can 
obtain the same result as for the constant smoothing length case. 

Figure [7] shows the dispersion relation calculated using the method of Appendix A in 
the one-dimensional case, where we use the equation of motion for variable smoothing 
length, Eq. (IMll . quintic spline interpolation, po = 0-8, po.eos = 10.0, and Pq = —9.2. The 
other parameters are the same as in Fig. [31 The results shown are for Csmooth = 1-0, 2.5, 
and 5.0. 
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Figure 7: Dispersion relation for variable smoothing length and negative pressure with po,eos = 10.0. Red 
squares show Csmooth = 1-0, green circles show Csmooth = 2.5, and blue triangles show Csmooth = 5-0. 
The solid curve shows the analytical dispersion relation for constant smoothing length. 


As shown in Fig. [71 for Csmooth = 1.0 the calculation becomes unstable at almost all 
wavelengths, but for Csmooth = 5.0 the calculation is stable at all wavelengths. Moreover, 
the solid line shows the analytical dispersion relation for constant smoothing length 
with the same parameters, and we confirm that this solid line almost coincides with 
the dispersion relation of variable smoothing length with Csmooth = 5.0. Therefore, 
we can obtain the same result as for constant smoothing length if we increase Csmooth- 
Appropriate values of Csmooth may depend on the equation of state, the other parameters, 
the spatial dimension, and so on. We should use sufficiently large Csmooth in practical 
calculations. 

We should note that even with the variable smoothing length, the tensile instabil¬ 
ity should not appear if we calculate the exact integration of convolution of EoM (e.g. 
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Eq. (77) of [ 13 ). Therefore, unwanted results for small Csmooth is due to the poor approx¬ 
imation of the convolutions with the variable smoothing length. We can qualitatively 
understand why the extension to the variable smoothing length itself aggravates the 
tensile instability in the negative pressure. The reason is as follows: if two particles ap¬ 
proach each other, the density becomes large and the smoothing length becomes short. 
This makes the shape of the kernel function sharp, and then the force between these 
two particles becomes strong. Therefore, if the pressure is negative, the attractive force 
becomes strong, and then this enhances the tensile instability. Large Csmooth makes this 
nature of the variable smoothing length weak for the perturbations of short wavelength. 

If Csmooth is larger, the number of neighbor particles required to obtain the correct 
smoothing length is larger. According to 17|, the number of neighbor particles for 
each dimension is bTyCsmooth for one dimension, for two dimensions, and 

for three dimensions. Especially for the two- and three-dimensional cases, 
the increase of computational cost with larger Csmooth is significant. However, from 
Eq. (1361) for the multidimensional case, the spatial variation of the smoothing length 
becomes smoother than the one-dimensional case, and then we can consider that the 
Csmooth required to obtain the same result as for constant smoothing length becomes 
smaller. 

This approach to variable smoothing length is not ideal. For example, if we improve 
the derivation of the equation of motion for variable smoothing length, we may formulate 
a method that works for Csmooth = 1- This may emphasize the importance of rigorous 
formulation of Godunov SPH for variable smoothing length. 


4-4- Summary of results 

Table □ shows a summary of stability for the number of dimensions (d), various 
interpolations, and positive and negative pressure in the constant smoothing length case. 



d= 1, P > 0 

d= 1, P < 0 

d = 2 or 3, P > 0 

d = 2 or 3, P < 0 

Linear 

0 

X 

0 

X 

Cubic 

0 

X 

X 

0 

Quintic 

X 

0 

X 

0 


Table 1; Summary of stability for the number of dimensions, various interpolations, and positive and 
negative pressure. A circle (O) shows stable interpolation and a cross (x) shows unstable interpolation. 
d represents the spatial dimension and P represents the pressure. 


According to the test calculation of the shock tube problem in [17| , the method with 
linear interpolation is not so accurate at the contact discontinuity. Therefore, it is good 
to use cubic spline interpolation for positive pressure and quintic spline interpolation 
for negative pressure in the one-dimensional case. However, for the two- and three- 
dimensional cases we must use linear interpolation for positive pressure. For negative 
pressure, both cubic spline interpolation and quintic spline interpolation are stable, but 
we recommend avoiding quintic spline interpolation because of its heavy computational 
cost. Thus, in the two- and three-dimensional cases it is good to use linear interpolation 
for positive pressure and cubic spline interpolation for negative pressure. 

To achieve conservation of all momentum and energy, we can use Pi + Pj as the 
criterion for the pressure sign. Therefore, the equation of motion in one dimension 
becomes 
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^.Lubic if(P,+P,)>0, 
^.Luintic if(P*+P,)< 0 , 


( 70 ) 


and for two and three dimensions 



(71) 


To achieve conservation of energy, we must use the same interpolation of for the 
equation of energy. In the case of variable smoothing length, we use Eq. (l34l) in the same 
way as Eq. and Eq. (ED, but we should use larger Csmooth to achieve the same results 
as in Table [H 

5. Test Calculation 

In this section, we perform test calculations to confirm that the results of the analyses 
in Section 4 are also valid in actual calculations. 

The time integration method is a simple predictor-corrector method as follows: first 
we calculate the acceleration at the nth time step, Vi^n, using the physical quantities at 
the nth time step, and then we calculate the time-centered position and velocity as 


At 


'^i,n+l/2 — '^i,n T 'di,n ^ ’ 



(72) 


Then we calculate the time-centered acceleration, Vi^n+l/ 2 ^ using the time-centered phys¬ 
ical quantities. Einally, we calculate the position and velocity at the next time step as 


— '^i,n U2,n+l/2^t, 

. 1 . .2 
n^n+l = -f Vi^nAt + -Vi^n+l/2At . 


(73) 


The time step At is determined by the Courant condition at each time step, 



( 74 ) 


In all the test calculations of this paper we use Ccfl =0.5. 
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5.1. Convergence test 

In this subsection, we perform a convergence test to check the accuracy of the Go¬ 
dunov SPH method. We calculate the sound wave in two dimensions as a problem for the 
convergence test. The density in the unperturbed state is set to 1.0. We use po.eos = 0.1 
and Cs = 1.0 for the equation of state. Thus the pressure in the unperturbed state is 
0.9. Simulations are performed in the square domain, x, y € [0.0,1.0], and we assume a 
periodic boundary condition. In the unperturbed state, positions of particles are 

given as a square lattice. The initial positions and velocities of particles are. 


Xi = Xi + (0.001/fc) sin(fca;i), 

Vi,x = —a;(0.001/fc) cos(fci^), 

Vi =¥i, 

Vi,v = 0, (75) 

where k = uj = 2 tt. In this case, the amplitude of the density perturbation is O.OOI. 
We consider the case of variable smoothing length with Csmooth = 1-0, and use linear 
interpolation for because only the linear interpolation is stable for case with the 
positive pressure in two dimensions. 

To measure the error, we calculate the difference of the density as. 


e 


1 


Ntot 

|Pref(*’i) 

i=l 




(76) 


where Wot is the total number of particles, Pref(^i) is a reference density at position r^. In 
this convergence test, we adopt the result with Wot = 512x512 as the reference density. 
We conduct the calculations with Wot = 16x16, 32x32, 64x64, 128x128, 256x256, and 
compare with the reference density. To reduce the error due to discrete time integration, 
time stepping At is set to very small value 5 x 10“^ for all resolution. We evaluate the 
error of density after 100 time-steps. In Fig.lH e is plotted as a function of the average 
smoothing length, which represents the resolution. 

As we can notice from Fig.[8j the error is proportional to h^. Therefore, it is confirmed 
that the Godunov SPH method has spatially second-order accuracy. 


5.2. Sound wave in one dimension 

First we calculate the propagation of a sound wave in the one-dimensional case as 
the simplest test problem. 

The particles in the unperturbed state are given equal spacing Ax = 0.01, and xi 
represents the unperturbed position. The initial positions and velocities of particles are 
given by 


Xi = Xi + O.OlAi sin(fca;i), 

Vi = —O.OlAxw cos(fc^), (77) 


where k is the wave number and oj is the angular frequency of the sound wave. In this 
calculation we use fc = 27r, and to set Cg = uj/k = 1 we use w = 27r. In other words, we 
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Figure 8; Result of convergence test of the Godunov SPH method for a sound wave. Horizontal axis 
shows the average smoothing length, and vertical axis shows the relative error e. Solid line represents 
the line oc h^. Note that the smoothing length is proportional to the average particle spacing. 


resolve one wavelength of the sound wave by 100 particles. We use Eq. (15^ for P*j and 
turn off the dissipation. 

All particles have the same mass m = 0.008, and then the average density is po = 0.8. 
We use the equation of state of Eq. (IMl) . In the positive pressure case we use po.eos = 0.6 
and the average pressure is Pq = 0 .2, and in the negative pressure case we use po,eos = 1-0 
and Pq = —0.2. The boundary condition is periodic. 

Figure[n]shows the density distribution of the sound wave with negative pressure using 
cubic spline interpolation and a constant smoothing length h = 0.01 immediately after the 
instability becomes noticeable {t = 1.154972). Figure [TO] shows the density distribution 
for the same conditions but using quintic spline interpolation at t = 5.249869. 

As we can see in Fig. [5] and Fig. [TUI our calculations remains stable in the case of quin¬ 
tic spline interpolation even at five times longer than the time at which the calculation 
with cubic spline interpolation becomes unstable. The calculation with quintic spline 
interpolation remains stable even at t = 50.0. 

In order to quantify the stability of these calculations, we calculate the time evolution 
of the minimum distance between all pairs of particles. We define this quantity as 

Tinin = min \ri-rj\/Ax. (78) 

all pairs of ij 

If the calculation is stable, rmin remains about unity because the distance to the nearest 
neighbor particle is about Ax. In contrast, if the particles begin to clump, rmin becomes 
small. Figure [11] shows the time evolution of r^nin in the case of linear interpolation, 
cubic spline interpolation, and quintic spline interpolation during the calculation of the 
negative pressure sound wave, and Fig.lTUlshows the same result but for positive pressure. 

As shown in Fig. [11] and Fig.[TU] the calculation for negative pressure with quintic 
spline interpolation remains stable, but for positive pressure the calculations with linear 
interpolation and cubic spline interpolation remain stable. These results are the same 
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Figure 9: Density distribution of the sound wave using cubic spline interpolation and a constant smooth¬ 
ing length for negative pressure at t = 1.154972. 
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Figure 10: Density distribution of the sound wave using quintic spline interpolation and a constant 
smoothing length for negative pressure at t = 5.249869. 


as those of Section 4. Here the lines of unstable calculations are broken because these 
calculations are terminated. 

Next, we consider the case of variable smoothing length. As in Section 4, we use 
Po,eos = 10.0, Pq = —9.2, and rj = 1.0. The other parameters and initial conditions are 
the same as in the case of constant smoothing length. Figure [T3] shows the time evolution 
of rmin when Csmooth = 1.0, 2.5, and 5.0. 

As shown in Fig. [131 in the calculations with Csmooth = 1-0 and Csmooth = 2.5 
the particles stick together soon after the beginning of the calculation. However, with 
Csmooth = 5-0 our calculation is stable. These results also agree with the linear stability 
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time 


Figure 11: Time evolution of rmin^ the minimum distance between adjacent particles, in the case of 
linear interpolation (red solid line), cubic spline interpolation (green dashed line), and quintic spline 
interpolation (blue dotted line) during the calculation of the sound wave for negative pressure. 



time 


Figure 12: Same as Fig. llll but for positive pressure. 


analyses. Therefore, stable calculation is possible if we make Csmooth sufficiently large in 
the actual calculation with variable smoothing lengths. 

Next, to show the steepening of the wave with high amplitude in the negative pressure 
case, we start from the initial conditions defined as 


Xi = Xi + O.lAx sin(fca;i), 

Vi = —0.1Aa;a;cos(A:^), (79) 
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Figure 13: Time evolution of rmin for the case of variable smoothing length with Cgmooth = 1-0 (red 
solid line), Cgmooth = 2.5 (green dashed line), and Cgmooth = 5.0 (blue dotted line). 


where all parameters are the same as for the case of Fig.[TUJ and we use quintic spline 
interpolation and a constant smoothing length. Figure [TJ] shows the density distribution 
at t = 0.0, t = 5.0, and t = 10.0. 
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Figure 14: Evolution of the density distribution of the sound wave with negative pressure. The red solid 
curve shows the density distribution at t = 0.0, the green dashed curve shows t = 5.0, and the blue 
dotted curve shows t = 10.0. 


As we can see from Fig.[T4l the wave profile shows gradual steepening as in the case of 
positive pressure. Therefore, even in a negative pressure medium, nonlinear steepening 
of a sound wave is expected to generate shock waves. 
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5.3. Shock tube problem 


The calculation of a wave is very simple and physical quantities do not vary largely. 
Thus, we should test our method for more nonlinear processes, such as those involving a 
discontinuity. In this subsection, we consider a shock tube problem in a negative pressure 
medium. To handle the shock wave we need small but finite physical viscosity. In the 
Godunov method, we use the Riemann solver to introduce a minimum but sufficient 
viscosity. Therefore, we use the analytical solution of the Riemann solver for Eq. (1551) 
that is derived in Section 3.3 to obtain P,* . 

I 

In [17l| , Inutsuka constructed a second-order Riemann solver by considering the spatial 
gradients of the physical quantities. We also use this second-order Riemann solver. The 
gradients of the physical quantites are calculated by Eq. (|S]). 

However, in the negative pressure region, the gradient of pressure for perturbation of 
Nyquist frequency can not be calculated correctly. Using Eq. (|5]), the gradients of density 
and pressure are calculated as. 


3 

vp. = Y. 


d 

'rnj-^Wir, - Vj.h), 

P d 

mj ——W{ri - rj,h). 
Pj dri 


(80) 

(81) 


Here, if we assume Nyquist frequency perturbation and equal mass particles, the 
density of all particles are the same, and the pressure of all particles are also the same in 
the case of equation of state Eq. (1551) . In that case, the gradients of density and pressure 
become, 


Vpi = —IU(ri - rj,h), 



rj,h), 


(82) 

(83) 


where, m, P and p are mass, pressure and density respectively. Thus, if P < 0, the 
gradients of pressure and density are anti-parallel. This is unphysical, because VP = 
CgVp with positive means that two gradients are parallel. 

If we use inconsistent gradient of pressure, we tend to estimate the resultant pressure 
of the Riemann problem for approaching particles mistakenly smaller, and this makes 
attractive force stronger. Therefore, we should calculate the gradient of pressure by a 
much better method. One possible way is. 


VP, = Cl,Vp,, (84) 

where, Cs,i is local sound speed of i-th particle, and Vpi is calculated by Eq. (ISO]) . In 
this paper, we use this modified gradient of pressure. In the case of variable smoothing 
length, we replace h of Eq. (IMl) with hi. 
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Here, we slightly modify the monotonicity constraint. We use a first-order Riemann 
solver when there are some particles that have opposite-sign gradients within their neigh¬ 
borhood. This condition is written for pairs of i-th and j-th particles as 


where 



(85) 


( 86 ) 


and / represents p or P. In this subsection, if there is any one particle j that satisfies 
the condition Eq. (1851) within 3/ii from the z-th particle, we use the first-order Riemann 
solver for the z-th particle. 

We use the equation of state given by Eq. (I!?!?l) . Cs = 1.0, and po.eos = 2.5. The initial 
discontinuity of the shock tube problem is at x = 0, and the initial parameters are 


Pl = 2.0, PR = 1.0, 

Pl = -0.5,Pr = -1.5, 

=0.0,z;,,R = 0.0, (87) 

where L denotes the physical quantities on the left-hand side of the initial discontinuity, 
and R denotes the right-hand side. The mass of each particle is m = 0.01, the particle 
spacing on the left-hand side is Axl = 0.005, and on the right-hand side it is Axr = 0.01. 
We put 200 particles on the left-hand side and 100 particles on the right-hand side. The 
boundary condition is a wall boundary {vx{x = — 1) = z;a;(x = 1) = 0). 

Eirst, we compare the results of cubic spline interpolation and quintic spline interpo¬ 
lation. Here we use a constant smoothing length h = 0.01. Eigure [15] shows the density 
distribution in the case of cubic spline interpolation immediately after the instability 
becomes visible {t = 0.4200), and Eig.lTBl shows the density distribution in the case of 
quintic spline interpolation at t = 0.4200. Here, the solid curve corresponds to the ana¬ 
lytical solution of the shock tube problem that is derived with the analytical solution of 
the Riemann problem in Section 3.3. 

As shown in Eig. ll5l and Eig. 1161 in the case of cubic spline interpolation the numerical 
instability occurs at the initial discontinuity, but in the case of quintic spline interpolation 
there is no instability. Note that the contact discontinuity does not exist because we use 
Eq. (1551) for the equation of state and the pressure only depends on the density. The 
result of this simulation matches the analytical solution well. 

Next, we test the calculation with variable smoothing length. Fig.[T7| shows the 
density distribution with Csmooth = 2.0 at t = 0.4200. Here we use rj = 1.0. 

Our test calculations show that the calculation with Csmooth = 1-0 becomes unstable, 
but the calculation with Csmooth = 2.0 is stable even with variable smoothing length. 

Therefore, the stability of the shock tube problem also agrees with the result of the 
linear stability analysis of a sound wave. 
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Figure 15: Density distribution of the shock tube problem with negative pressure at t = 0.4200 obtained 
by cubic spline interpolation. The solid line corresponds to the analytical solution. 



Figure 16: Density distribution of the shock tube problem with negative pressure at t = 0.4200 obtained 
by quintic spline interpolation. The solid line corresponds to the analytical solution. 


5.4- Sound wave in two dimensions 

In this subsection, we consider a sound wave in the two-dimensional case. In the 
unperturbed state, the particles are put on a square lattice with a side length of Aa; = 
0.01. The initial positions and velocities are 
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Figure 17: Density distribution at t = 0.4200 obtained by quintic spline interpolation and variable 
smoothing length with Csmooth = 2.0. 


Xi= Xi+ 0.01 Aic sm{kxi), 

Vi^x = —O.OlAiw cos(fc^), 

Vi = ¥i, 

Vi,y = 0 . ( 88 ) 

In this problem, we consider a sound wave that propagates toward the x-direction. 
We use k = 27r/(40Ax) = S.Ott and w = S.Ott to set Cs = w/fc = 1.0. Thus we resolve 
1 wavelength with 40 particles. The mass of each particle is m = 0.8 x 10““^ and the 
average density is po = 0.8. The equation of state is again given by Eq. (l3^ . and in the 
positive pressure case we use po,eos = 0.6 and Pq = 0.2, and in the negative pressure case 
we use po.eos = 1-0 and Pq = —0.2. We use Eq. (15^ for Pt. The boundary condition is 
periodic for both the x- and y-directions. In this calculation we use a constant smoothing 
length h = 0.01. 

Figure [TS] shows the density distribution for the x-direction wave with negative pres¬ 
sure using cubic spline interpolation at t = 3.004534. In the two-dimensional case with 
negative pressure, the calculation using linear interpolation is unstable. However, as we 
can see in Fig. 1181 the calculation with cubic spline interpolation remains stable, and we 
confirm that this remains stable at least until t = 50.0. 

Figure [12] shows the evolution of Tmin during the calculation with negative pressure 
using linear interpolation, cubic spline interpolation, and quintic spline interpolation, 
and Fig. 1201 shows with positive pressure. We can see that rmin becomes small in the 
calculation with linear interpolation for negative pressure, and cubic spline interpolation 
and quintic spline interpolation for positive pressure. The result for cubic spline inter¬ 
polation is in contrast to the one-dimensional case. These results also agree with the 
results of the linear stability analyses of a sound wave. 
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Figure 18: Density distribution in the i-direction for two-dimensional sound wave propagation with 
negative pressure at t = 3.004534 using cubic spline interpolation. 



Figure 19: Time evolution of rmin during the calculation of two-dimensional sound wave propagation 
with negative pressure.The red solid line shows linear interpolation, the green dashed line shows cubic 
spline interpolation, and the blue dotted line shows quintic spline interpolation. 


The Cartesian lattice is not an ideal configuration compared to more stable config¬ 
uration like the densest-sphere packing. Thus we conduct the same test for the case of 
a regular triangular lattice. All conditions other than the position in the unperturbed 
state are the same as the previous test. For simplicity, we test linear interpolation and 
cubic spline interpolation, and we do not test quintic spline interpolation, because the 
stability of cubic and quintic spline interpolation is the same in two dimensions. Figure 
M shows the evolution of rmin using linear interpolation and cubic spline interpolation. 

As shown in Fig. 1211 the calculations with linear interpolation for positive pressure and 
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Figure 20: Same as Fig. llOl but for positive pressure. 



Figure 21: Time evolution of rmin during the calculation of two-dimensional sound wave propagation. 
In the unperturbed state, particles are put on a regular triangular lattice. The red solid line shows linear 
interpolation with positive pressure, the green dashed line shows cubic spline interpolation with positive 
pressure, the blue dotted line shows linear interpolation with negative pressure, and the black chain line 
shows cubic spline interpolation with negative pressure. 


cubic spline interpolation for negative pressure are stable, but the calculations with linear 
interpolation for negative pressure and cubic spline interpolation for positive pressure are 
unstable. These results are consistent with the case of the Cartesian lattice. 

Moreover, we conduct the calculation of the sound wave with different direction of 
the wave number vector. In the unperturbed state, the particles are put on a Cartesian 
lattice. We use fe = (S.Ott, S.Ott) and w = \k\ = 5.0-\/27r. The initial positions and the 
velocities are. 
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ri = rl + Sr^, 

Sn = (O.OlAa;, O.OlAa;) sin(fc • ri), 

Vi = —a;(0.01Aa;, O.OlAa;) cos(fc • ri). (89) 

The other parameters and conditions are the same as in the case with the wave 
number vector is along the x-direction. Figure [H] shows the evolution of rmin in this case 
using linear interpolation and cubic spline interpolation. 



Figure 22: Time evolution of rmin during the calculation of two-dimensional sound wave propagation 
with k cc (1,1). The red solid line shows linear interpolation with positive pressure, the green dashed line 
shows cubic spline interpolation with positive pressure, the blue dotted line shows linear interpolation 
with negative pressure, and the black chain line shows cubic spline interpolation with negative pressure. 


As we can notice from Fig. [211 the results are the same with those of the lattice-grid 
aligned cases. 

5.5. Shock tube problem in two dimensions 

In this subsection, we consider a shock tube problem in two dimensions. To demon¬ 
strate the validity of our method in a situation where the sign of pressure changes spatially 
and temporally, we calculate a shock tube problem with different signs of pressure in left 
and right side of initial discontinuity. 

Also in this subsection, we use the equation of state given by Eg. (1551) . Cs = 1.0 
and po.eos = 2.5. We set the initial discontinuity at x = 0, and the region for x < 0 
corresponds to the left side of the initial discontinuity, x > 0 corresponds to the right- 
hand side. The initial parameters are, 


Pl = 4.0, PR = 1.0, 

Pl = 1.5,Pr = -1.5, 

Vx,-L = I’y.L = O.OjUx.R = Uy,R = 0.0. 
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We put SPH particles on a regular triangular lattice. The side length of this regular 
triangle is, 0.01 for the left side and 0.02 for the right side. The mass of each particle is 
m = 2-\/3 X 10“^. In this subsection, we use constant smoothing length with h = 0.02 
for simplicity. The boundary condition for the cc-direction is a wall boundary {vx{x = 
— 1) = Vx{x = 1) = 0), and for the y-direction we use a periodic boundary. In the same 
way as in one dimension, we use the second-order Riemann solver for elastic equation of 
state. 

Figure [53] shows the density distribution at t = 0.6000 when we use appropriate 
interpolation of as Eq. (ED, and Fig. 1331 shows the same density distribution at the 
same time but we use only linear interpolation independent of the sign of pressure. 



Figure 23: Density distribution of the shock tube problem in two dimensions at t = 0.6000 obtained by 
choosing appropriate interpolation. The solid line corresponds to the analytical solution. 


As shown in Fig. [531 if we use appropriate interpolation depending on the sign of 
pressure, we can calculate without any problem. However, as shown in Fig. 1241 if we 
ignore the sign of pressure and we use only linear interpolation, the particles at contact 
surface make clustering and we can not calculate correctly. Therefore, our method is 
valid even in the case where the sign of pressure changes spatially and temporary. 

5.6. Equilibrium test with random noise 

In this subsection, we calculate the time evolution from a noisy initial condition in 
three dimensions. This test is very similar to that of Q. Initially, the particles are 
put on a face-centered cubic lattice with nearest-neighbor distance dnn = 0-1-^. We 
add random noise of the position with the amplitude of O.Idnn to each direction of each 
particle. Initial velocities of all particles are set to 0. We put 4000 particles in the 
computational domain, and assume a periodic boundary condition. 

We set the mass of each particle m = 2.5 x 10“^. Thus the density in unperturbed 
state is po = I-O- We use equation of state Eq. (l33ll . and Cg = I-O, po.eos = 0.8 for positive 
pressure, po.eos = 1-2 for negative pressure. We assume constant smoothing length with 
h = dnn, and we use Riemann solver for Eq. (1331) to achieve the equilibrium state. 
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Figure 24: Same as Fig. l23l but we use only linear interpolation. 


Following the test of Q, we calculate the time evolution of q^m 


which is defined as, 


^min — mill \Vi ^jl/dnn- 

all pairs of ij 


(91) 


This ^min is an indicator for the regularity of particle distribution. If the particles are 
put on a perfect face-centered cubic lattice, gmin becomes one, while if the particles are 
clustered this value is close to 0. In a typical grass-like distribution (a uniform-density 
equilibrium distribution), gmin ~ 0.7. Figure ESl shows the time evolution of gmin in 
the case of linear interpolation and cubic spline interpolation for positive and negative 
pressure respectively. 

As shown in Fig.[25j gmin of linear interpolation for negative pressure and cubic spline 
interpolation for positive pressure become close to 0 during the time evolution, which 
means that the particles are clustering. On the other hand, gmin of linear interpolation for 
positive pressure remains almost initial value. Moreover, that of cubic spline interpolation 
for negative pressure increases gradually and becomes close to I, which means that the 
distribution of particles gradually settles onto a face-centered cubic lattice. Therefore, 
these results in three-dimensional simulation also agree with those of linear stability 
analysis. 


6. Summary and Future Work 

In the SPH method, there is a numerical instability that results in the clustering 
of particles under certain conditions, and this tensile instability is significant in the 
negative pressure regime. In this paper we show that the formalism of the Godunov SPH 
method [l7| is viable for mitigating the tensile instability. We formulate higher-order 
approximations for the interpolation of V (s) and conduct linear stability analyses for the 
various equations of motion of the Godunov SPH method. 


36 






1 


0.8 

. . 

0.6 


C 


E 

\ \ 

O' „ , 

\ 

0.4 

H 

0.2 

K 




0.1 1 ' 10 

time 


Figure 25: Time evolution of ^min for th® simulations starting from noisy initial conditions in three 
dimensions. The red solid line shows linear interpolation with positive pressure, the green dashed line 
shows cubic spline interpolation with positive pressure, the blue dotted line shows linear interpolation 
with negative pressure, and the pink chain line shows cubic spline interpolation with negative pressure. 


We conclude that in the one-dimensional case the stable interpolations are linear 
interpolation and cubic spline interpolation for positive pressure media, and quintic spline 
interpolation for negative pressure media. In the two- and three-dimensional cases, linear 
interpolation is stable for positive pressure, and cubic spline interpolation and quintic 
spline interpolation are stable for negative pressure. In the case of variable smoothing 
length, calculations with sufficiently large Csmooth remain stable. Therefore, we can 
suppress the tensile instability by an appropriate order of interpolation without any 
additional and artificial terms. Additionally, we have derived the analytical solution of 
the Riemann problem for the equation of state given by Eq. (I33p , and we confirmed that 
our analytical solution agrees with the result of the numerical simulation of the shock 
tube problem. 

In practical calculations of elastic dynamics, we need to formulate how to handle 
deviatoric stress tensor that corresponds to the non-diagonal parts of the stress tensor 
(e.g., [13) 0|)' The extension of the present method to the non-diagonal stress tensor will 
be studied in our next paper. 
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Appendix A 

In Appendix A, we explain how to obtain the dispersion relation from the acceleration 
term. We assume that the wave number vector is along the x-direction. 

First, we define the unperturbed positions of SPH particles as 

r\ = = {lAx,mAx,nAx), i = l,2,---N (Al) 

where N shows the number of particles, and I, m, n denote integers 1, 2, •••. In 
other words, SPH particles are put on a square lattice with a side length of Ax in the 
unperturbed state. 

We consider the perturbed positions to be 


n = {xi +dxi,yi,Zi), 

Sxi = e exp[i(A:^ — ujt )], (A2) 

where k is the wave number of the perturbation, ui is the angular frequency of the 
perturbation, e is the amplitude, and i not in subscript shows the imaginary unit. 
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We expect that the acceleration can be expressed as 


^ ^ (5xi,0, 0). (-^^) 

In practice, we can calculate by taking the ratio of ax,i and Sxi'. 

=-ax,il5xi. (A4) 

Thus we can obtain the dispersion relation by taking the ratio of the displacement 
and the acceleration for various wave numbers k: 


U! 


1 ^ 

2=1 


^x,i \ 

Sxi) 


(A5) 


Appendix B 

In Appendix B, we describe the linear stability analysis of the Godunov SPH method. 

We assume particle positions as in Eq. dsa. Here we consider e as an inhnitesimal 
constant and neglect second or higher orders of e. Only the x-component of the specific- 
volume gradient and the acceleration appear, and the other components vanish because 
the perturbation is only along the x-axis. We assume infinitely accurate time integration, 
and ignore the effect of the discretization of time integration. 

We assume the masses of all particles, m, are the same. We use the equation of motion 
given by Eq. dSOl), and use Eq. ([^ for T’* . For , we use Eq. ([Ml) for linear interpolation, 
Eq. for cubic spline interpolation, and Eq. for quintic spline interpolation. 

First, we linearize the density of particle i using Eq. ([5]): 


_ ^ ^ ^ ^ 

p, = '^mW{ri - rj,h) Ri mW {r\ - r], h) + '^m{Sxi - Sxj)-^W{j\ -r],h) 


r Q 1 

= 'y^ mW{ri — r], h) + — exp[—ik{3^—3^)])-^zzW{ri— rj,h) Sxi- (Bl) 


The first term on the right-hand side of Eq. (ED is the density of the unperturbed state. 
This term is almost the same as the average density po- The terms with odd functions 
ofxi—H^ vanish when we sum over subscript j. Thus, Eq. (IBII) becomes 


Pi ~ Po(l - iD6xi), (B2) 

where D is defined by Eq. (EZl). From Eq. (IB2I) . we can immediately find the linearized 
specihc volume as 


Vi =—{l + iDSxi). (B3) 

Po 

Next, we linearize the ^-component of the gradient of the specific volume. From 
Eq. (1551) . the gradient of the specific volume becomes 
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dxi pI dxi 

1 Q2 1 

~-? m{Sxi — Sxj)-—^W{ri — rj,h) = - CpSxi, (B4) 

Pq j OXi Po 

where Cp is defined by Eq. (1571) . Finally, we linearize the second-order derivative of the 
specific volume using Eq. (I55|) : 


d'^V, 

d4 




——CpDbxi- 

Pq 


The linearized pressure of particle i is 


(B5) 


P^ = Po + 5P^Po + Cl5p = Po- iClpaD5x,. (B6) 

We substitute these linearized physical quantities into Eq. (1501) . and we define the 
coefficients a, b, and so on using Eq. (ED. Finally, we obtain the analytical solution of 
such as Eq. Eq. ((Ml) , and Eq. (IHHl) . 
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